import numpy as np
import pandas as pd
from sklearn.metrics import r2_score  #评价回归预测模型的性能
import matplotlib.pyplot as plt #画图
import seaborn as sns
import time

df = pd.read_csv("FE_boston.csv")
print(df.head())

##从原始数据分离 特征X 和输出y
X = df.drop(['MEDV','log_MEDV'], axis=1)
y = df['MEDV']
feat_names = X.columns
print(feat_names)

'''
当数据量比较大时，可用train_test_split从训练集中分出一部分做校验集；
 样本数目较少时，建议用交叉验证。
 在线性回归中，留一交叉验证有简便计算方式。
下面将训练数据分割成训练集和测试集，
只是让大家对模型的训练误差、校验集上的测试误差估计、和测试集上的测试误差做个比较。
'''
#讲数据分割训练数据和测试数据
from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,random_state=33,test_size=0.2)
print(X_train.shape)

##线性回归
from sklearn.linear_model import LinearRegression
#1.使用默认配置初始化学习器
lr = LinearRegression()
#2.用训练数据训练模型参数
lr.fit(X_train, y_train)
#3.用训练好的模型对测试集预测
y_test_pred_lr = lr.predict(X_test)
y_train_pred_lr = lr.predict(X_train)

# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性
fs = pd.DataFrame({"columns":list(feat_names), "coef":list((lr.coef_.T))})
fs.sort_values(by=['coef'],ascending=False)
print(fs)

# 使用r2_score评价模型在测试集和训练集上的性能，并输出评估结果
#测试集
print('The r2 score of LinearRegression on test is', r2_score(y_test, y_test_pred_lr))
#训练集
print ('The r2 score of LinearRegression on train is', r2_score(y_train, y_train_pred_lr))

###在训练集上观察预测残差的分布，看是否符合模型假设：噪声为0均值的高斯噪声
f, ax = plt.subplots(figsize=(7,5))
f.tight_layout()
ax.hist(y_train-y_train_pred_lr,bins=40,label='residuals Linear', color= 'b', alpha=.5)
ax.set_title('histogram of residuals')
ax.legend(loc='best')
plt.show()
##观察预测值和真值的散点图
plt.figure(figsize=(4, 3))
plt.scatter(y_train, y_train_pred_lr)
plt.plot([-3, 3], [-3, 3], '--k')   #数据已经标准化，3倍标准差即可
plt.axis('tight')
plt.xlabel('True price')
plt.ylabel('Predicted price')
plt.tight_layout()
plt.show()

##
# 正则化的线性回归,(L2正则：岭回归)
from sklearn.linear_model import RidgeCV
print("岭回归:")
#1.设置超参数
alphas = [0.01, 0.1, 1, 10,100]
#2.创建ridgeCV实例
ridge = RidgeCV(alphas=alphas, store_cv_values=True)
#3 训练模型
ridge.fit(X_train, y_train)
# 预测
y_test_pred_ridge = ridge.predict(X_test)
y_train_pred_ridge = ridge.predict(X_train)
#评估
print('r2 score of rigde on test:',r2_score(y_test, y_test_pred_ridge))
print('r2 score of ridge on train:', r2_score(y_train, y_train_pred_ridge))

#可视化
mse_mean = np.mean(ridge.cv_values_, axis = 0)
plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas),1))

#这是为了标出最佳参数的位置，不是必须
plt.plot(np.log10(ridge.alpha_)*np.ones(3), [0.28, 0.29, 0.30])

plt.xlabel('log(alpha)')
plt.ylabel('mse')
plt.show()
print ('alpha is:', ridge.alpha_)
##查看各个特征的权重系数, 系数的绝对值大小反应特征的重要性
fs = pd.DataFrame({'columns':list(feat_names),
                   'coef_lr':list(lr.coef_.T),"coef_ridge":list(ridge.coef_.T)})
fs.sort_values(by=['coef_lr'], ascending=False)

print('特征系数：',fs)


'''
线性回归：L1正则->Lasso
'''
##Lasso/ L1正则
# class sklearn.linear_model.LassoCV(eps=0.001, n_alphas=100, alphas=None, fit_intercept=True,
#                                    normalize=False, precompute=’auto’, max_iter=1000,
#                                    tol=0.0001, copy_X=True, cv=None, verbose=False, n_jobs=1,
#                                    positive=False, random_state=None, selection=’cyclic’)
from sklearn.linear_model import  LassoCV
print('Lasso/L1正则:')
#1. 设置超参数搜索范围
#alphas = [ 0.01, 0.1, 1, 10,100]

#2. 生成学习器实例
#lasso = LassoCV(alphas=alphas)

#1. 设置超参数搜索范围
#Lasso可以自动确定最大的alpha，所以另一种设置alpha的方式是设置最小的alpha值（eps） 和 超参数的数目（n_alphas），
#然后LassoCV对最小值和最大值之间在log域上均匀取值n_alphas个
# np.logspace(np.log10(alpha_max * eps), np.log10(alpha_max),num=n_alphas)[::-1]

#2 生成LassoCV实例（默认超参数搜索范围）
lasso = LassoCV()

#3 训练
lasso.fit(X_train, y_train)

#4.测试
y_test_pred_lasso = lasso.predict(X_test)
y_train_pred_lasso = lasso.predict(X_train)
# 评估，使用r2_score评价模型在测试集和训练集上的性能
print('r2 socre of lassoCV on test is ', r2_score(y_test, y_test_pred_lasso))
print('r2 socre of lassoCV on train is ', r2_score(y_train, y_train_pred_lasso))
##可视化
mses = np.mean(lasso.mse_path_, axis=1)
plt.plot(np.log10(lasso.alphas_), mses)
# plt.plot(np.log10(lasso.alphas_)*np.ones(3), [0.3, 0.4, 1.0])
plt.xlabel('log(alpha)')
plt.ylabel('mse')
plt.show()

print('alpha is:', lasso.alpha_)
# 看看各特征的权重系数，系数的绝对值大小可视为该特征的重要性
fs = pd.DataFrame({"columns":list(feat_names), "coef_lr":list((lr.coef_.T)), "coef_ridge":list((ridge.coef_.T)), "coef_lasso":list((lasso.coef_.T))})
fs.sort_values(by=['coef_lr'],ascending=False)
print('lasso 特征系数：\n', fs)
